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1. Introduction 

The main goal of particle physics today is the search for new physics. A primary focus of that 
search is the precision determination of the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements 
or, more accurately, an overdetermination. To achieve this, many flavor-changing processes are 
studied by both experiment and theory and the results are combined to look for inconsistencies, 
which could indicate new physics at work [1, 2, 3]. At the same time, as hints of new physics 
appear; beyond the Standard Model phenomenologists use the information to constrain, support, or 
rule out new models [4, 5, 6]. In both cases, lattice-QCD calculations of hadronic matrix elements 
are crucial. 

Neutral B mixing is a promising place to look for new physics for two reasons. First, in the 
Standard Model it occurs via a loop process in which the contributions are suppressed both by 
CKM matrix elements, and, for all but the top quark, masses of the loop quark. So, any new 
physics could present a relatively large signal. Second, many recent reviews of the persistent 2-3a 
tension between the Standard Model and flavor physics experiments indicate that new physics in B 
mixing is a likely explanation [1, 2, 3, 4, 5]. 

The most general AB = 2 effective hamiltonian contains eight dimension-six operators 

5 3 

^eff = £Q^; + £C,-^M (1.1) 

!=1 (=1 

with 

ff2 = {b''Lq") (PLql^), = {b'^Lq^) {PRq""), (1.2) 

^3 = (^«L/) {b^Lq''), 

where a and j8 are color indices, L and R are the projection operators ^(1 it 75), and 1^1,2,3 are 
obtained from (^1.2.3 by L — > /?. 

Parity conservation of QCD ensures the matrix elements for ^1^2,3 are physically equivalent to 
those for ^1,2,3, leaving five independent matrix elements. The matrix element i^i appears 
in the Standard Model expression for the B meson mass difference AM^. The matrix element 
(5^1^21^^) is needed for renormalization of (Bjj|^i|B^), and all three can be useful in studies of 
the width difference AF and its constraints on new physics [7, 8]. The matrix elements (B°|i^4|B°) 
and (B^|(^5|B°) are needed for generic extensions beyond the Standard Model (some examples are 
given in Refs. [6]). 

In the Standard Model, 

(^1^) %(Ai)IW/(B^l^i|fi^>(M), (1-3) 

where the quantities in parentheses and tjb are known factors, Vij is a CKM matrix element, and we 
use relativistic normalization for the states. Measurements of Mdq have sub-percent eiTors [9, 10], 
so our ability to constrain the CKM matrix contribution |Vi/,Vf* | is limited by how precisely we 
know (B^l^ijB^). Current, published lattice-QCD calculations [11] report errors of just under 
14% on (B^|^i|B°). 
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For mixing only, an unquenched calculation of {B^\ 1^1,2,3 1^^) appeared in Ref. [12], though 
it includes data at only one lattice spacing, a 0.12 fm. The only calculation to date of ^4|B^) 
and (B^l^sjB^) appeared a decade ago [13]. In that work all five matrix elements were calculated 
using the quenched approximation for the sea, static heavy valence quarks with an interpolation to 
ntb, and a linear extrapolation in the light valence mass to reach Bj. These approximations are no 
longer required in modern lattice-QCD calculations. We include three sea quarks to simulate up, 
down, and strange; work directly at m^; include data at two or more lattice spacings; and use rooted, 
staggered chiral perturbation theory to guide a simultaneous chiral-continuum extrapolation. This 
work presents the first full error budget for the beyond the Standard Model matrix elements. 

Historically, it has been useful to parametrize the matrix elements as [7, 14] 

(sO|^,-|S«)(At) = Q Ml fl B'ili^), (1.4) 

with coefficients c,- = (2/3, ^Vi2, 1/12, 1/2, i/e). The B^^, known as bag parameters, character- 
ize deviation from the so-called vacuum saturation-approximation [15]. Today, it is possible to 
calculate the {B^^\^i\B^^) without appealing to this approximation and, therefore, without need for 
introducing the bag pai^ameters. Nevertheless, because of their historical use, at a later stage in 
this project, we intend to combine results for with our collaboration's results for fs^ 

to extract the bag parameters in a manner that accounts for coiTclations between the two calcula- 
tions. In this work we present preliminary results for the matrix elements and refer the reader to 
Refs. [16, 17], and references therein, for values of /b^ calculated with similar fermion actions and 
gauge configurations. 

In Sec. 2 we describe details of our calculation. A discussion of the full error budget appears 
in Sec. 3. In Sec. 4 we report preliminary results, and in Sec. 5 we summarize and discuss the 
remaining work to be done. 

2. Calculation 

In this section we review the generation of coiTelation function data, the extraction of matrix 
elements from fits to the data, and the renormalization and chiral-continuum extrapolation of the 
matrix elements. 

2.1 Generating Correlation Function Data 

Correlation function data are generated using the MILC gauge configurations [18] with 2 + 1 
flavors of asqtad staggered [19] sea quarks. Bottom valence quarks are simulated using the Fer- 
milab interpretation of the Sheikholeslami-Wohlert (clover) action [20] with hopping parameter Kh 
tuned to produce the observed B^ meson mass [21]. Light valence quarks are simulated with the 
asqtad action. We work in the meson rest frame. 

Table 1 lists the ensembles used to date and selected valence quark parameters. For each of the 
"conf configurations we average over ?isrc sources with temporal spacing T /4 and spatial distribution 
randomized to reduce coiTclations. At this stage of the calculation, different subsets of data have 
been analyzed for the different operators. The subsets used are noted in Table 1. 
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0.0982 


1996 X 4 






X 


0.06 
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Table 1: Summary of the MILC ensembles and valence masses used in this analysis. Each lattice has 
an approximate spacing a, spatial extent L, and temporal extent T . The tadpole improvement factor is 
determined from the expectation value of the average plaquette, mq = (plaque tte)'/^. The light and strange 
sea-quark masses are w/ and w/,, respectively. The sets of valence masses used, niq, generally run from 
m/,/10 to m/„ with values given by valo.12 = {0.005, 0.007, 0.01, 0.02, 0.03, 0.0349, 0.0415, 0.05}, valo.09 
= {0.0031, 0.0047, 0.0062, 0.0093, 0.0124, 0.0261, 0.031}, and valo.oe = {0.0018, 0.0025, 0.0036, 0.0054, 
0.0072, 0.016, 0.0188}. The bottom quark hopping parameter is fQ,. An "X" marks ensembles used in the 
the analysis of each operator 

Data are generated as ensemble average estimates of the vacuum expectation values 

c2p'(0 = £(b°(^x)B«(0,0)A (2.1) 

and 

C-^\h.t2) = £ (b°(?2, X2) 0) B°(?i, Xl)). (2.2) 

Xl,X2 

The B meson creation operator is given by 

x)^ =Y.W{t, y)5(x,y) 75 q{t, x), (2.3) 
y 

and the four-quark mixing operator is 

^,•(0,0) = v/(o, 0) rp'^(o, 0) v/(o, o) rf ^^(o, o). (2.4) 

Propagators arising from the Wick contraction of naive light-quark fields q are related to staggered 
quark propagators, generated in our simulation, following Ref. [22]. The heavy-quark field i// is a 
four-component Dirac spinor, smeared in the B meson creation operator to improve overlap with 
the ground state. The smearing function 5(x,y) is based on the quarkonium IS wavefunction [23]. 
Heavy-quark fields in the mixing operator are rotated to remove leading order discretization er- 
rors [20]. The T\ ' represent the spin structures appearing in the quark bilinears of Eq. (1.2). The 
three-point correlation function setup is illustrated in Fig. 1 . 
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Figure 1: A is created at rest at ti < 0. At time 0, via mixing operator it oscillates into a B^, which 
is subsequently annihilated at t2 > 0. The heavy (light) quark propagator is indicated by a double (single) 
Une. To extract the hadronic mixing matrix element we need the three-point correlation function depicted 
here and the two-point correlation functions corresponding to the propagation of the B^ from ti to and of 
the B° from to f2. 



Although parity conservation of QCD ensures the hadronic mixing matrix elements of ^1,2,3 
are physically equivalent to those of ^1,2,3, the finite statistics estimates of Eq. (2.2) need not 
generate data that are numerically identical. We average these physically equivalent data to increase 
statistics associated with the determination of the matrix elements of ^1,2,3- 

2.2 Fitting Correlation Function Data 

We model two- and three -point coiTclation function data by projecting the vacuum expectation 
values of Eqs. (2.1) and (2.2) onto bases of and B° energy eigenstates, 

C^P\t)= £ + (2.5) 

cfP'(?i,?2) =""1'— f^J^la^MB^] (_i)'"('.+i)+«fe+i)^-£„^2^-£„,^,^ (2.6) 



n,m=Q 



where oscillating, opposite-parity state contributions [22] are incorporated. The effects of periodic 
boundary conditions are negligible in fitted three-point data and are omitted from Eq. (2.6). 

We extract the ground state quantity {Bq\ffi\Bq) /Mb^^ by simultaneously fitting the two- and 
three-point correlation function data to Eqs. (2.5) and (2.6). In principle the sums run over infinitely 
many radial excitations of the pseudoscalar B meson, in practice we must fit with finite N'^'^^K We 
also restrict data included in the fit (fmin < ? < fmax)- Fits using data at short times must account for 
excited state contributions by including adequate numbers of states. Despite added difficulty, the 
relatively clean signal in the data at short times makes it desirable to include them in the fit. The use 
of a Bayesian fitting routine [24, 25] allows us to deal with the difficulty of fitting numerous states 
and incorporate data at short times. Increasing f^ax has the benefit of utilizing more data. However, 
this also introduces an increasing level of noise and can lead to a poorly determined covariance 
matrix. 

From plots of correlation function data, scaled to remove leading exponential decay, we de- 
termine a range of potential f^^'" and t^ax^ values to study. ^ The range of potential fmax" values 
is determined by the onset of significant noise. For we consider from t^^^^ = 2 until excited 



'For C^^^ (f 1 , f2 ) , we could choose distinct and f2.i 
symmetry of data generated by Eq. (2.2) prompts the simplifications: = fi.min = h.min and /^^x = fi.niax = tl.n 



5 



Er^ and B^. mixing 



CM. Bouchard and E.D. Freeland 



state contributions have significantly decreased. We perform fits witli A/^^'^p' = 2,4, and (some- 
times) 6? Starting with the two-point correlation functions, we fit each combination of A/^^p', t^^, 
and fmax and select representative fits at each N^^^ from a common plateau of fit results for Zq and 
£0 versus t^^^ and fmax> thereby ensuring stability with respect to N^p\ t^^, and f^ax- Then, fixing 
^min ^^'^ ^max Separately for each A/^^p', we repeat the procedure for a simultaneous fit of the two- and 
three -point correlation function data. For these fits, we fix N^^^ = N^^^ and determine the optimum 
N^^\ ^m'in' ai^d ^max- This Systematic procedure, outlined in more detail in [26, 27], allows us to 
obtain fit results that are stable in the number of states and data included in the fit. We find the fit 
results are largely insensitive to changes in fmax'"- 

2.3 Renormalization 

Lattice results for the hadronic mixing matrix elements are combined with perturbatively cal- 
culated Wilson coefficients, i.e. C,- and Q of Eq. (1.1), to obtain physical quantities. Therefore, the 
lattice and continuum results must be brought to the same scheme and scale. We use the mass of 
the bottom quark to set the scale and match to the continuum MS— NDR scheme. At one loop the 
matrix elements mix under renormalization, 

{B\ ^,-|5)^-N™(m,) = £ [5ij + a.,{q*K,-j{amb)] {B\ ffj\B)'^\ (2.7) 

where the matching coefficients give the difference between the one-loop lattice and continuum 
renormalizations [28]. We evaluate as{q*), as in [29], in the "V scheme" [30] at the scale q* = 2 /a. 
MS-NDR is not sufficient to define the matrix elements of operators ^2 and ^3 in dimensional reg- 
ularization; one must also choose a set of "evanescent" operators. The renormalization coefficients 
and must be calculated accordingly. Below, we present results for two schemes, which we 
label "BBGLN" [31] and "BJU" [32]. The one-loop conversion between the two schemes was first 
given in Ref. [13]. A more general transformation of the two-loop anomalous dimension matrix 
can be found in [33]. For convenience, we define the shorthand (^,) = {B\&i\B)^^^^^^{mi,), for 
the matrix elements after extrapolating the renormalized results to the continuum. 

2.4 Chiral-Continuum Extrapolation 

The continuum, next-to-leading order (NLO), chiral perturbation theory for hadronic mixing 
matrix elements [34], modified [35] for the rooted, staggered, light-quark actions (rS;^PT), guides 
a simultaneous extrapolation to physical quark mass and the continuum. The rS;^PT expression is 
fit to the renormalized matrix elements of Eq. (2.7) to determine the low energy constants. Effects 
from heavy-light meson Havor- and hyperfine-splittings [16] and finite volume (expected to be 
< 1%) are not explicitly incorporated in the chiral extrapolations presented here, but will be in the 
final stage of analysis. Chiral fits are performed using a Bayesian fitting routine [24, 25] with priors 
constraining each fit parameter. Additional details will be given in future publications. 

A chiral expansion including only NLO terms fits the data reasonably well when data included 
in the fit are limited to light valence masses {i.e., values of lUq such that the pseudoscalar meson 

^We use equal numbers of parity-even and parity-odd states in our fits. 
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mass of Eq. (2.9) satisfies (M^^ ri)^ < 1).^ To include data at all simulated valence masses we find 
it necessary to include NNLO analytic terms. There are seven potential NNLO analytic terms [36] 



aX.' a\2Mf,+Ml), M^^^{2Mfi+Mf 



m: 



{2Mf,+Ml)\ 



a\ {2Mti+Mtn), (2.8) 



where the (squared) pseudoscalar meson mass is given by 



Mf:=Bo{mi+mj) 



(2.9) 



and Bq is a low-energy constant from ;tPT. Each of these terms comes with a coefficient to be 
determined by the fit. At a minimum, we find it necessary to include terms proportional to a~M^^ 
and M^^{2Mfj +Mf^jj) in order to achieve suitable fits to data at all valence masses. 

Reported central values are obtained from those fits including all seven NNLO terms. The 
simultaneous chiral-continuum extrapolations of (^,) /Mb are shown in Figs. 2, 3, 4, 5, and 6. The 
p- values of these fits are very near one. This is not the case for NLO fits or fits including various 
subsets of the possible NNLO analytic terms of Eq. (2.8). Inflated p- values, at least partially, result 
from the fact that we use the Bayesian augmented and degrees of freedom. In addition, the size 
of our statistical sample may be limiting our ability to resolve small eigenvalues of the coiTelation 
matrix. We are studying this issue. 



X7dof = 10.0/50 p-value =1.0 
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Figure 2: {^i): chiral-continuum extrapolation data and fits. (Red/orange) squares are 0.12 fm data, (green) 
circles are 0.09 fm data, and (blue) x 's are 0.06 fm data. The extrapolation/interpolation is the black line 
with a shaded (turquoise) error band. Black bursts mark the point of extrapolation (interpolation) for the B^/ 
(Bs) meson. 



'For reference, at the strange valence mass, {M^q ri)^ ~ 1.3 
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Figure 3: {^i}- chiral-continuum extrapolation data and fits in the BBGLN scheme described in Sec. 2.3. 
(Red/orange) squares are 0.12 fm data, (green) circles are 0.09 fm data, and (blue) x's are 0.06 fm data. 
The extrapolation/interpolation is the black line with a shaded (turquoise) error band. Black bursts mark the 
point of extrapolation (interpolation) for the B^j (B^) meson. 
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Figure 4: (^3): chiral-continuum extrapolation data and fits in the BBGLN scheme described in Sec. 2.3. 
(Red/orange) squares are 0.12 fm data and (green) circles are 0.09 fm data. The extrapolation/interpolation 
is the black line with a shaded (turquoise) error band. Black bursts mark the point of extrapolation (interpo- 
lation) for the B^j (By) meson. 
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0.8 I ^ ^ ^ 1 

0.0 0.5 1.0 1.5 2.0 

Figure 5: (ftj): chiral-continuum extrapolation data and fits. (Red/orange) squares are 0.12 fm data, (green) 
triangles are 0.09 fm data. The extrapolation/interpolation is the black line with a shaded (turquoise) error 
band. Black bursts mark the point of extrapolation (interpolation) for the B,/ (Bs) meson. 




0.0 0.5 1.0 1.5 2.0 




Figure 6: (i^^'s): chiral-continuum extrapolation data and fits. (Red/orange) squares are 0.12 fm data, (green) 
triangles are 0.09 fm data. The extrapolation/interpolation is the black line with a shaded (turquoise) error 
band. Black bursts mark the point of extrapolation (interpolation) for the B^i {B^) meson. 
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Source of Error [%] (iffi) 








„ statistical 8.6 
chiral-continuum systematic 12 


6.8 16 
11 3.3 


4.3 
0.2 


5.5 
4.4 


statistical 6.7 
chiral-continuum systematic 1.8 


4.6 10 
6.6 4.5 


2.5 
1.6 


3.4 
3.7 



Table 2: Statistical and chiral-continuum systematic errors for the hadronic mixing matrix elements are 
determined from the data by methods outlined in Sec. 3. Here we show results from the use of the BBGLN 
evanescent operators in the renormalization, as described in Sec. 2.3. To the precision reported here, errors 
in the BJU scheme differ only for the chiral-continuum systematic error for mixing via {i?2), where we 
found the error to be 6.4%. 



Sourceof Error [%] (^,) 

scale (n) 3 

Kb tuning 4 

light-quark masses 1 
heavy-quark discretization 4 

one-loop matching 8 

finite-volume effects 1 

subtotal 10 



Table 3: Estimates of errors based on our collaboration's previous works [16, 37]. We anticipate significant 
reduction in the leading sources of error once all ensembles at a w 0.06 fm and at a w 0.45 fm have been 
fully incorporated into the analysis. 

3. Errors 

The error budget given in this work is complete in that it accounts for all known, significant 
sources of error in our lattice-QCD calculation. We have split the list of en^ors into two tables. 
Table 2 lists the statistical and chiral-continuum eiTors, which ai^e determined directly from the 
analysis presented here, for each (^,) individually. EiTors estimated from our collaboration's pre- 
vious works [16, 37] are listed in Table 3 and are taken (for now) to be the same for all (^,). EiTors 
associated with omitting charm and heavier sea-quark contributions, isospin breaking, and elec- 
tromagnetic effects are expected to be insignificant in comparison to the errors in Tables 2 and 3. 
Further study, and an ongoing analysis of a 0.06 fm and a 0.045 fm data, should lead to a 
reduction of our leading sources of error in Tables 2 and 3. 

Statistical errors are obtained by propagating a distribution of 600 bootstrap ensembles of 
con^elation function data {i.e., Eqs. (2.5 and 2.6)) through the con^elation function and rS;t;PT fitting 
procedures, providing non-parametric error estimates that account for congelations among the data. 
This en^or also includes uncertainty associated with chiral fit parameters /;;: and gB*B7i by setting 
prior widths based on parameter uncertainties. 
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The systematic error associated with the chiral-continuum series truncation and extrapolation 
is obtained by varying the rS;^PT fit ansatz. For both the S° and matrix elements, the addition 
of the two minimum required NNLO analytic terms of Eq. (2.8) is necessary, and sufficient, to 
achieve good fits to all valence mass data that are stable against the addition of remaining NNLO 
terms. We estimate the systematic eiTor to be the shift in fit central values between fits with the two 
required NNLO analytic terms and those with all NNLO analytic terms. 

The remaining errors, quoted in Table 3, are conservative estimates based on previous col- 
laboration publications. The errors associated with scale-setting and jq, tuning are estimated from 
Ref. [16]. Entries for light-quark mass, heavy-quark discretization, and one-loop matching errors 
are estimated from Ref. [37]. Finite volume error estimates are based on Ref. [37] and inflated to 
accommodate an expected increase [38] from hyperfine- and flavor-splitting effects. 

4. Results 

Results from our chiral-continuum extrapolations are in ri units. We convert to GeV using the 
lattice determination r\ = 0.31 17(22) fm [16] and the conversion factor 1 GeV fm = 5.0677312(1). 
In Ueu of the matrix elements, in Table 4 we report the phenomenologically more useful quan- 

7 (i) 

tity /g^ B^' . Results ai^e given in the continuum MS— NDR scheme at scale mi, using both the 
BBGLN [31] and BJU [32] choices of evanescent operators. Applying a heavy-quark relation 
between the matrix elements of ^i, 1^2, and ^3, the conversion [13] for bag parameters can be 
expressed as 

Bg B,j 4jj. 1^5^? ^9 / ' ^ ^ 

^i3)BJU ^ j^i3)BBGLN /_4^(1) + ^j^iA . (4.2) 

For the results given in Table 4, we bring the matrix elements to the BBGLN or BJU scheme before 
the chiral-continuum extrapolation. Converting a posteriori with Eqs. (4.1) and (4.2) yields values 
in very good agreement. 

The scheme dependence of the /^^ B^^ from Table 4 matches the corresponding scheme de- 
pendence in the Wilson coefficients with which they will be used. In particular, one must use the 
BBGLN or BJU values for (^2) and (1^3), according to the choice of evanescent operator scheme 
used in the Wilson coefficient calculation. As mentioned in Sec. 2.3, (^1), {^4), and ((^5) ai^e 
independendent of this choice. 

Errors are obtained by adding in quadrature the appropriate errors from Table 2 and the subtotal 
of Table 3. As shown in Table 2 and illustrated in Fig. 4, the matrix element (B^|^3|Bg) suffers 
from much larger errors than the others and is more difficult to analyze. We are optimistic that 
the addition of data at finer lattice spacings, and refinements in analysis techniques, will lead to 
improvements. 

5. Summary and Outlook 

We report on our ongoing, 2+1 flavor, lattice-QCD calculation of the hadronic contribution 
to S*' and mixing. Table 4 gives preliminary results for the hadronic mixing matrix elements for 
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[GeV^] 


Bd 

BBGLN BJU 


s 

BBGLN BJU 


J Bg Bq 


041 UTS') 






0.0574(92) 0.0538(87) 


0.086(11) 0.080(10) 


-' Bq Bq 


0.058(11) 0.058(11) 


0.084(13) 0.084(13) 


fl 

■' Bq Bq 


0.093(10) 


0.135(15) 


fl 

■' tSq Bq 


0.127(15) 


0.178(20) 



Table 4: Preliminary results for f\B^g^ — (^,) / CiMj^, evaluated in the continuum MS— NDR scheme at 
scale wifo. "BBGLN" and "BJU" denote the choice of evanescent operators in the renormalization schemes 
of Refs. [3L 32], respectively. For / = 1,4,5, the schemes are the same. Although the chiral-continuum 

9 (3) 

extrapolation was done separately for both schemes, its output for fsq^s^ is the same for the number of 
digits shown. The meson masses, needed in combination with the quantity extracted in Eq. (2.6), are 
known to < 0.01% [10] and do not contribute significantly to the error. 

a basis of four-quark, dimension-six, AB = 2 mixing operators that spans the space of all possible 
hadronic contributions in the Standard Model and beyond. At the current stage of analysis our 
errors are competitive with published Standard Model matrix element results. For beyond the 
Standard Model matrix elements, this is the first unquenched calculation and the first update in ten 
years. 

Additional data at finer lattice spacings, a « 0.06 fm and a « 0.045 fm, is currently being 
analyzed. We anticipate that its inclusion, and a more rigorous analysis of the errors of Table 3, will 
result in a significant reduction in the error budget. The discretization and one-loop matching eiTors 
of Table 3 ai^e based on a 0.09 fm data, and will decrease when re-evaluated for a w 0.045 fm. 
Recent data runs will allow for a new analysis of K"/, tuning errors as well. We are finalizing the 
calculation of the one-loop perturbative matching coefficients of Sec. 2.3 and plan to incorporate 
finite volume effects and hyperfine- and flavor-splittings in the chiral-continuum extrapolation. 
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